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I .   Introduction 

This  report  describes  position  estimation  procedures  that 
are  based  on  models  which  relate  positional  uncertainty  to 
various  measurement  and  estimation  errors.   The  procedures  were 
developed  to  be  used  in  analyses  that  relate  the  effect  of  posi- 
tional uncertainty  on  tactical  performance  through  such  factors 
as  weapon  accuracy  as  well  as  to  be  used  operationally. 

In  the  models,  positions  are  on  a  plane  surface  (flat  earth 
model).   Because  of  this  condition,  the  models  are  not  intended 
for  use  in  situations  where  the  earth's  figure  is  significant. 
In  addition,  positional  errors  are  determined  by  independent 
normally  distributed  random  variables  with  known  means,  variances 
and  covariances.   The  support  for  this  condition,  other  than  its 
mathematical  convenience,  is  that  it  has  been  used  by  others, 
for  example,  see  Reference  1.   To  use  the  models,  one  is  required 
to  specify  the  means,  variances  and  covariances  of  the  error 
random  variables.   For  these  models,  this  can  be  done  by  speci- 
fying the  systematic  errors  (biases)  and  the  error  magnitudes 
(standard  deviations) . 

The  procedure  that  is  described  in  Section  II  relates  position 
estimates  that  are  based  on  bearings  on  or  from  stations  to 
station  position  uncertainty.   One  application  for  the  model  is 
the  analysis  of  the  effect  of  sonobuoy  position  uncertainty  on 
position  estimates  determined  with  passive  directional  sonobuoys . 

The  procedure  that  is  described  in  Section  III  relates  posi- 
tion estimates  that  are  based  on  lines  of  position  to  line  of 


position  uncertainty.   In  addition  to  estimating  the  uncertainty 
in  position  estimates  based  on  celestial  observations,  the  pro- 
cedure could  be  used  to  determine  error  ellipses  for  LORAN  fixes 
if  the  standard  deviation  values  required  by  the  procedure  could 
be  obtained. 

The  procedure  that  is  described  in  Section  IV  can  be  used 
to  combine  position  estimates  from  various  sources.   The  proce- 
dure which  is  based  on  conditions  that  should  not  be  too  restric- 
tive  in  most  cases  provides  both  a  composite  position  estimate 
and  error  ellipse. 


II.   Station  Position  Uncertainty  and  Position  Estimates 

The  procedure  that  is  described  in  this   section  relates 
position  estimates  that  are  based  on  bearings  on  or  from  stations 
to  station   position  uncertainty.   The  procedure  is  based  on  a 
model  that  is  an  extension  of  one  that  is  described  in  Appendix 
1.   The  model  is  defined  as  follows:   Each  station  position  error 
is  determined  by  an  independent  bivariate  normal  distribution 
with  a  zero  mean  vector  and  a  known  covariance  matrix.   Observed 
bearing  lines  on  or  from  a  station  are  parallel  to  true  bearing 
lines.   The  distance  of  each  observed  bearing  line  from  its  true 
bearing  line  is  determined  by  an  independent  normal  distribution 
with  a  zero  mean  and  a  standard  deviation  a. 

Because  a  station's  position  error  is  determined  by  a  bivar- 
iate normal  distribution,  the  perpendicular  distance  s  between  a 
line  at  the  assumed  location  of  an  observed  bearing  line  and  the 
observed  bearing  line  is  determined  by  a  normal  distribution. 
The  relation  between  the  bivariate  normal  distribution  that 
describes  the  station's  position  uncertainty  is  indicated  in 
Figure  1.   In  the  figure,  the  positive  y-axis  direction  is  north, 
the  positive  x-axis  direction  is  east,  and  the  origin  of  the 
coordinate  system  is  at  the  assumed  station  position.   The  x'y'- 
coordinate  system  is  oriented  so  that  the  positive  x'axis  is 
coincident  with  the  major  axes  of  the  elliptical  contours  of 
the  bivariate  normal  distribution  that  determines  the  station 
position  error  and  so  that  the  bearing  6  of  the  positive  x'-axis 


assumed  location  of 
the  observed  bearing  line 


Figure  1.   The  geometry  associated  with  the  determination  of  os, 
the  standard  deviation  of  the  normal  distribution 
that  determines  the  distance  between  the  assumed 
location  of  an  observed  bearing  line  and  the  observed 
bearing  line.   The  assumed  (mean)  position  of  the 
station  is  at  the  origin.   The  ellipse  represents  a 
contour  on  the  probability  density  surface  of  the 
bivariate  normal  distribution  that  determines  the 
station's  position. 


satisfies  the  condition:  0°    <$  <  180°.   The  angle  6  is  an 
observed  bearing.   The  x"y "-coordinate  system  is  oriented  so 
that  the  positive  y"-axis  is  in  the  direction  of  the  observed 
bearing  line.   As  a  consequence  of  these  relationships, 

2       2     2  2     2 

og   =  ax,   sin  (6-9)  +  a  ,   cos  (6-e) 

2         2 
where  a  ,   and  a  ,   are  the  elements  of  the  station  position 
x        y 

error  covariance  matrix  relative  to  the  x ' y ' -coordinate  system. 

A  procedure  for  determining  position  estimates  that  are 
based  on  bearings  on  or  from  stations  at  known  positions  is 
described  in  Appendix  1  of  this  report.   The  procedure  is  based 
on  a  model  that  is  equivalent  to  one  that  is  described  in  Ref- 
erence 2.   The  procedure  that  is  described  in  this  section  is 
based  on  a  model  that  is  an  extension  of  it. 

The  procedure  in  Appendix  1  of  this  report  is  based  on  a 
model  in  which  a  station's  bearing  error  is  determined  by  a 
normal  distribution  with  zero  mean  (bias)  and  standard  devia- 
tion e.   The  bearing  error  is  related  to  the  distance  on  a  cir- 
cular arc  between  a  station's  true  bearing  line  and  observed 
bearing  line.   The  arc  is  on  the  circle  with  its  center  at  the 
station  that  passes  through  an  initial  estimate  of  an  object's 
position.   This  distance  is  determined  by  a  normal  distribution 
with  mean  zero  and  standard  deviation  a  =  re  where  r  is  the 
range  of  the  initial  estimate  from  the  station  and  the  standard 
deviation  e  is  measured  in  radians.   In  the  model,  arc  distance 


is  approximated  using  a  first  order  approximation  which  in 
effect  replaces  the  circle  with  its  tangent  line  at  the  initial 
estimate.   As  a  consequence,  the  distance  u  on  the  tangent  line 
between  the  observed  bearing  line  and  the  true  bearing  line  is 
determined  by  a  normal  random  variable  with  mean  zero  and  stan- 
dard deviation  o.   The  distance  u  can  be  expressed  in  terms  of 
w,  the  distance  on  the  tangent  line  between  the  observed  bearing 
line  and  the  initial  estimate  that  is  also  determined  by  a  normal 
random  variable  with  standard  deviation  a,  and  v,  the  distance 
on  the  tangent  line  between  the  true  bearing  line  and  the  initial 
estimate.   And,  as  shown  in  Appendix  1  of  this  report,  this 
distance  can  be  expressed  in  terms  of  the  unknown  coordinates 
of  the  object's  position. 

The  effect  of  station  position  uncertainty  is  accounted  for 
by  the  distance  s  between  the  observed  bearing  line  and  the 
assumed  observed  bearing  line.   In  the  model  s  is  determined  by 

a  normal  distribution  with  mean  zero  and  standard  deviation  a 

s 

as  given  above.   This  approximation  is  consistent  with  the  first 
order  approximation  of  arc  distance.   As  a  consequence  of  these 
two  approximations,  all  of  the  bearing  lines  are  replaced  by 
lines  parallel  to  the  line  joining  the  initial  estimate's 
position  and  the  station's  assumed  position,  both  of  which  are 
known  positions.   The  geometry  involved  in  shown  in  Figure  2. 
The  modified  relationships  resulting  from  the  introduction  of 
station  position  uncertainty  are  shown  in  Figure  3.   The 
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Figure  2.   The  geometry  of  the  position  estimation  model.   The 
bearing  (3  and  range  r  of  the  initial  estimate  from 
the  assumed  station  position  have  the  role  of  3 
and  r  in  Appendix  1  of  this  report  • 


object 


assumed  location 
of  the  observed  bearing  line 


observed  bearing 
line 


w'  =  w  -  s 


Figure  3.   The  quantities  u,  v  and  w  correspond  to  the  quantities 
u,  v  and  w  in  Appendix  1  of  Reference  3.   The  auxiliary 
quantities  u1  and  w1  are  defined  in  the  figure.   The 
replacement  of  a^  by  a^  +  os^    in  the  procedure  in 
Appendix  1  of  this  report  is  justified  by  noting  that 
with  this  substitution  u1  and  w*  are  equivalent  to  u 
and  w. 


important  thing  to  note  is  that  u1  =  u  -  s  is  determined  by  a 

normal  distribution  with  mean  zero  and  standard  deviation 

2      2    h 
[o      +  o   ]   and  that  otherwise  U'  is  equivalent  to  U  with 

respect  to  the  procedure  in  Appendix  1  of  this  report.   As  a 

consequence  of  this,  the  procedure  can  be  extended  to  include 

2      2  h 
station  uncertainty  by  replacing  a  by  [a   +  a   ]2  where  ever 

it  is  used.   In  this  case,  a  =  re  where  r  is  the  range  of  an 

initial  estimate  of  an  object's  position  from  the  assumed  station 

position  and  e  is  the  bearing  error  (standard  deviation)  in 

radians  of  the  bearings  associated  with  the  station. 


III.   Line  of  Position  Uncertainty  and  Position  Estimates 

The  procedure  that  is  described  in  this  section  relates 
position  estimates  that  are  based  on  lines  of  position  to  line 
of  position  uncertainty.   The  procedure  is  based  on  a  model  that 
is  defined  as  follows:   Lines  of  position  are  straight  lines. 
Observed  lines  of  position  are  parallel  to  true  lines  of  posi- 
tion.  The  distance  of  an  observed  line  of  position  from  a  true 
line  of  position  is  determined  by  an  independent  normally  dis- 
tributed random  variable   with  known  mean  and  standard  deviation 
Lines  of  position  are  specified  in  terms  of  a  rectangular  coor- 
dinate system  with  the  origin  at  a  reference  point  as  shown  in 
Figure  4.   For  celestial  navigation,  an  appropriate  choice  for 
the  reference  point  would  be  the  assumed  position.   Since  bear- 
ing lines  are  lines  of  position,  this  model  differs  from  the 
model  that  is  described  in  Section  II  only  in  terminology. 
However,  operationally  the  use  of  the  model  that  is  described 
in  this  section  differs  in  the  way  the  standard  deviation  of  the 
distance  of  the  line  from  a  true  line  is  determined.   The  stan- 
dard deviation  associated  with  each  line  of  position  must  be 
specified.   If  this  can  be  done,  the  procedure  can  be  used. 
As  an  example,  suppose  that  the  values  are  a,  for  the  first 
line  of  position  and  a-  for  the  second  where  a,  >  c^    and 
that  the  lines  of  position  intersect  at  a  90°  angle.   In  this 
example,  the  minimum  area  confidence  (probability)  region  is 
an  ellipse  that  is  centered  on  the  estimated  position.   And  the 
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position 


initial  estimate 
and  final  estimate 


1st  line  of  position 


2nd  line  of  positic 


Figure  4.   The  geometry  when  only  two  lines  of  position  are  used 
with  the  procedure  in  Appendix  1  of  this  report.   The  lines  of 
position  correspond  to  bearing  lines  for  observed  bearings  9i  and 
9?  with  respect  to  the  procedure.   The  procedure  determines   the 
size  and  the  orientation  of  elliptical  confidence  regions  of  con- 
fidence p,  the  location  of  the  initial  estimate  and  the  location 
of  the  final  estimate.   The  locations  of  the  estimates  correspond 
when  only  two  lines  are  used.   If  the  lines  were  obtained  from 
sextant  observations,  the  assumed  position  should  be  the  reference 
position   in  which  case  the  lines  of  position  would  be  determined 
by  the  azimuth  anges  a,  and  a~  and  the  distances  p,  and  p~. 


11 


estimated  position   is  at  the  intersection  of  the  bearing  lines. 

2 
For  a  confidence  (probability)  of  containment  of  1  -  exp(-k  /2) , 

the  major  axis  of  the  ellipse  is  coincident  with  the  first  bear- 
ing line  and  it  is  of  length  2ko, ,  the  minor  axis  of  the  ellipse 

is  coincident  with  the  second  bearing  line  and  it  is  of  length 

2 
2ko2  and  the  area  of  the  ellipse  is  Trk  cj,a?  . 


12 


IV.   A  Composite  Position  Estimate 

The  procedure  that  is  described  in  this  section  is  for 
combining  position  estimates  for  an  object  from  independent 
sources.   It  is  based  on  the  following  model:   The  rectangular 
coordinates  of  each  position  estimate  are  determined  by  an  inde- 
pendent bivariate  normal  distribution  whose  covariance  matrix 
is  known  but  whose  mean  vector  is  not  known.   The  components  of 
the  mean  vector  for  each  of  the  distributions  are  the  unknown 
coordinates  x  and  y  of  the  object.   This  model  implies  that  the 
natural  logarithm  of  the  likelihood  function  for  a  set  of  n 
estimates  can  be  expressed  as  follows: 

n  ,   -1 

log  L  =  K  -  1/2  Z   (x.  -  x)  Z-  (x.  -  x) 

i=l  _1       L 

where  K  is  a  constant,  x.  is  an  estimate  vector  with  components  x- 

and  y.,  x  is  the  common  mean  vector  with  components  x  and  y, 

the  unknown  coordinates  of  the  object,  and  Z.  is  the  covariance 

2      2 
matrix  with  elements  a . G    ,  a.~   and  a.~~  .   The  maximum  likeli- 

l  yi  x  y        lxy 

hood  estimates  x  and  y  of  the  unknown  coordinates  x  and  y  are 
the  solutions  of  the  two  simultaneous  linear  equations  deter- 
mined by 


d  (log  L) 
33E 


x=x 

y=y 


=  o 


and 


9 (log  L) 

3y 


x=x 

y=y 


=  o 
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The  equations  can  be  written  as 

A  x  +  B  y  =  D 
B  x  +  C  y  =  E 

and  their  solution  as 

CD    -    BE                                     AE    -    BD 
x    =    7  y    =    y 

AC    -    B  AC    -    B 

where   A=    la.,    B   =    lb.,    C   =    Ec .  ,    D   =    Z(a.x.    +    b.y.), 
1  1  1  11  ri    ' 

E   =    Z(b.x.    +   c.y.),    a.    =    a.*    /d .  ,    b.    =   -    a.~-/d.,    c.    =    a • «    /d .  , 
i    l  iJ i  i  ly        i        i  lxy      11  ix   '     i' 

2         2  2 

d.  =  a.£  o.-  -  a.--   and  all  of  the  sums  are  for  i  from  1  to  n. 

Since  they  are  linear  combinations  of  the  estimates  x.  and 

y • ,  the  estimates  x  and  y  are  determined  by  a  bivariate  normal 

distribution.   Consequently,  all  that  is  required  to  determine 

this  distribution  is  its  mean  vector  with  components  y-  and  y- 

x      y 

2     2 

and  its  covariance  matrix  with  elements  a~  ,  a~   and  o~~.   The 

x  '   y       xy 

mean  vector  is  determined  by 

y-  =  E((CD  -  BE) /(AC  -  B2)}  =  x 
and 

UJ  =  E{(AE  -  BD)/(AC  -  B2)}  =  y 

And,  the  covariance  matrix  is  determined  by 

a-2  =  E  { (CD  -  BE)  -  E  (CD  -  BE) }2/(AC  -  B2)2 

=  (C2(F+I  +  2L)  -  2CB(H+K+M)  +  B2 (G+J+2N)  }/ (AC  -  B  )2 
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a-2  =  £{(AE  -  BD)  -  E(AE  -  BD)}2/(AC  -  B2)2 

=  {A2(G+J+2N)  -  2AB(H+K+M)  +  B2 (F+I+2L) }/ (AC  -  B2)2 

a--  =  E{[ (CD-BE) -E (CD-BE) ]•[ (AE-BD) -E (AE-BD) ] }/ (AC  -  B2)2 
xy 

=  { (AC+B  ) (H+K+M)-CB(F+I+2L)-BA(G+J+2N) }/ (AC  -  B2)2 

2    2  2    2  2  2    2 

where  F  =  Za.  a.-   ,  G  =  Zb.  a.-  ,  H  =  Za.b.a.~   ,  I  =  Zb.  a.-  , 
l   ix  l   ix  l  l  ix  i   ly 

2    2  2  2 

J  =  Zc .  a.-  ,  K  =  Zb.c.a.^  ,  L  =  Za.b.a.--,  M  =  Z(a.c.+b.  )a.~~ 
i   ly         l  l  ly         11  lxy         l  l   i    lxy 

and  N  =  Zb.c.a.~~  where  all  of  the  sums  are  for  i  =  1  to  n . 
l  l  lxy 

By  using  arguments  given  in  Appendix  1  of  this  report,  one 
can  show  that  the  axes  of  the  elliptical  confidence  regions 
associated  with  x  and  y  are  coincident  with  an  x' y '-coordinate 
system  where  the  transformation  from  the  xy-coordinate  system 

to   this  system  is  the  coordinate  axes  rotation  through  the 

2      2 
angle  y  defined  by  tan  2y  =  2a~V(a~   -  a~  )  .   For  a  confidence 
=    ■  J  '  xy    y      x 

p,  from  Reference  3,  the  minimum  area  confidence  region  is  an 

2 
ellipse  with  semi-axes  k  o~  ,  and  k  o-   ,  and  area  irk  o~,o~, 

x         y  i  ^   y 

where  k  =  [-2  ln(l-p)] 2  , 

2      2     2  .         2.2 

a~,   =  a-   cos  Y  -  2a~~  cos  y  sin  y  +  a~   sin   y, 
x'      x  xy  y 

2      2     2  2     2 

a~.   =  a~   sin  y  +  2a~~  cos  y  sin  y  +  a~   cos   y 
y1      x  xy  y 

and  the  center  of  the  ellipse  is  at  the  point  (x,y).   In  this 
coordinate  system,  °~t,}i  =  0  • 
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The  above  equations  can  be  used  to  specify  a  composite  posi- 
tion estimate  in  terms  of  the  location,  orientation  and  size  of 
an  elliptical  confidence  region  which  is  generally  the  form  in 

which  position  estimates  of  the  kind  that  are  being  considered 

2      2 
here  are  specified.   But,  since  values  of  a~   ,  a~  and  a^ 
t-  x     y       xy 

are  required  for  each  of  the  n  estimates  that  is  being  combined, 
a  way  is  needed  for  determining  these  values  given  the  orienta- 
tion and  size  of  an  elliptical  confidence  region.   A  procedure 
to  do  this  when  the  orientation  is  given  in  terms  of  the  direc- 
tion 6  of  the  major  axis  and  the  size  is  given  in  terms  of  the 
lengths  SMJ  and  SMI  of  the  semi-major  and  semi-minor  axes  and 
the  confidence  p  is  described  next. 

By  using  an  xy-coordinate  system  in  which  the  positive 
y-axis  direction  is  north  and  the  positive  x-axis  direction  is 
east  and  with  the  convention  0°  <  6  <  180°,  the  dependence  of 
the  value  of  the  rotation  angle  Y  and  of  the  order  relation 

between  a~ ,  and  a~ ,  on  the  value  of  the  major  axis  direction  6 
x       y 

is  indicated  by  the  following  table: 


0°  <  6  <   45°:     Y  =  6         and  a-,  >  a- , 

45°    6  <  135°:     y   =  6  -   90°  and  o~,  >  o~, 

x'     y ' 

135°  <  6  <  180°:     y   =  6  -  180°  and  a-,  >  a~, 

With  an  order  relation  and  a  value  for  p,  values  for  a~ ,  and 
a* ,    can  be  determined  with  values  for  SMJ  and  SMI.   With  values 

y 

2      2 

for  a~,  ,  a~.  and  y,    values  for  a~   ,  a^   and  a~~  can  be  found 
x'     y  x     y       xy 

from  the  following  equations: 
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2       2     2         2     2 

o~  =  a~    cos   y  +  a~    sin   y, 

x     x '  y  •        '  ' 

2       2     2         2     2 

a~  =  o~    sin   y  +  0%   cos   y 

y     x1  y 


and 


2       2 

a--  =  (o~    -  a~   )  sin  y  cos  y 
xy     y      x1 

2 
which  can  be  obtained  by  inverting  the  equations  above  for  a- , 

2 

and  o- , 

y 

As  an  example,  suppose  the  data  in  the  following  table 
represent  three  independent  position  estimates: 


X 

y 

6 

SMJ 

SMI 

k 

1st 

-3.7 

18.1 

59° 

36 

20 

2 

2nd 

11.8 

8.4 

105° 

37 

11 

2 

3rd 

0 

0 

146° 

45 

23 

2 

Here,  distances  are  in  nautical  miles  and  p  =  .86  in  each  case 
For  this  example  with  values  in  square  nautical  miles: 

a,-2  =  264.53,  a,-2  =  159.42  and  a   --  =  98.89 
lx  ly  lxy 

o„~2  =  321.35,  a  -2  =   51.15  and  a  ~-  =  -78 
2y  2y  2xy 

a.-2  =  249.20,  a  ~2  =  389.30  and  a,--  =  -173.4 
3x  '   3y  3xy 

These  values  give  the  following  composite  estimate: 
x  =  -2.46  nautical  miles  and  y  =  12.26  nautical  miles. 
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For  this  case,  a-  =  10.28  nautical  miles,  a~  =  6.21  nautical 
miles  and  o~~   =  -55.1  square  nautical  miles.   And,  for  k  =  2: 
SMJ  =  20.56  nautical  miles,  SMI  =  12.42  nautical  miles  and 
6  =  119°.   The  composite  confidence  region  and  its  three  com- 
ponent confidence  regions  are  shown  in  Figure  5. 

As  a  second  example,  suppose  each  position  estimate  is  deter- 
mined by  a  circular  normal  distribution.   Then  a.~  =  o.~  =  a. 
1  ix     ly     l 

and  cj.££  =  0  for  i  =  1  to  n.   In  this  case,  the  composite 

ixy 

estimate  is: 

x  =  (Zxi/o.)/(Zl/oi)  ,  y  =  (Zyi/ai)/(Zl/oi)  , 

2  2      2  2 

a~   =  n/(El/a.)   ,  o~   =  n/(£l/a.)   and  a~~  =  0. 
x  l    '   y     '    '  l        xy 

In  this  example,  since  x  and  y  are  determined  by  a  circular 
normal  distribution,  the  minimum  area  confidence  regions  are 
circles  and  orientation  is  not  an  issue. 
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y  North 


East 
x 


Figure  5.   The  ellipses  define  the  confidence  regions  of 
the  first  example.   The  composite  confidence  region  is  in 
black.   The  position  estimates  are  at  the  center  of  the 
ellipses.   The  numbers  indicate  the  order  of  the  estimate 
in  the  table  on  Page  17. 
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APPENDIX  1:   A  Position  Estimation  Procedure 

Consider  a  target  whose  position  is  unknown  and  a  set  of 
reference  stations  whose  positions  are  known.   Assume  conditions 
are  such  that  observed  target  bearings  from  the  stations  can  be 
considered  to  be  values  of  independent  normal  random  variables 
whose  means  are  equal  to  the  true  target  bearings  and  unknown  but 
whose  standard  deviations  are  known.   The  position  estimation 
procedure  that  is  described  here  is  a  maximum  likelihood  estimation 
procedure  that  is  based  on  this  assumption  and  the  assumption  that 
the  conditions  are  such  that  the  stations  and  the  target  can  be 
considered  to  be  located  on  a  plane  tangent  to  the  earth's 
surface  at  a  point  in  their  neighborhood. 

With  the  above  assumptions,  the  likelihood  of  observed 

bearings  6,,  9~,  ...,  8   from  stations  labeled  1,  2,  ...,  n  is: 
12         n 

Me,,  e,,  ...,  e  )  =  n —   exP  -  \    z  (e  -  $.)2/e2 

12         nl/2?e.  '      1    * 

l 

where  d> ,  ,  <fu ,  ...,  <J>   are  the  unknown  station  true  bearings  and 
1    2         n 

en,  e„,  ....  e   are  the  known  stations  standard  deviations. 
1    2         n 

To  a  first  order  approximation,  the  set  of  bearing  estimates 

<{>,.  du,  ...,  d)   that  are  determined  by  the  procedure  make 
12         n 

L(9  ,  9  ,  ...,  0  )  a  maximum  subject  to  the  constraint  that  the 
bearing  lines  determined  by  a  set  of  bearing  estimates  must  all 
pass  through  a  common  point.   The  common  point  is  the  estimate 
of  the  target's  position. 
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North 
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Initial 
Estimate 


Target 


Observed 
Bearing  Line 


Reference 
Station 


x 


last 


Reference  Point 

Figure  6.  The  coordinate  geometry.  The  coordinates  of  the  initial 
estimate  are  (x*,y*).  In  the  development,  the  reference 
point  is  at  the  initial  estimate. 
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In  order  to  impose  the  constraint  on  the  estimate  $ ,  through  |  , 

consider  the  quantities  u.  =  r.(0.  -  <!>.).  v.  =  r .  (d> .     -    3  .  )  and 

i     11     11     11     i 

w.  =  ri(°::  "  3.)  where  i  represents  any  station  number  from  1  to  n. 

Here,  <j) .     is  the  target's  true  bearing,  3.  is  the  bearing  of  an 

initial  estimate  of  the  target's  position,  0.  is  an  observed  target 

bearing  and  r.  is  the  range  of  the  initial  estimate  from  the 

station.   It  is  also  the  radius  of  a  circle  that  is  centered  on 

the  station  and  passes  through  the  initial  estimate  as  shown  in 

Figure  6  and  u.,  v.  and  w.  are  related  arc  lengths  on  this  circle 
3  11       i 

with  u.  =  w.  -  v. .   In  this  relation,  w.  is  known  and  v.  can  be 
ill  l  l 

expressed  in  terms  of  the  unknown  coordinates  of  the  target  with 

an  approximation  that  does  not  involve  <p .  .      To  do  this,  consider 

a  rectangular  coordinate  system  whose  origin  is  at  the  position 

of  the  initial  estimate  and  whose  axes  are  oriented  like  those 

shown  in  Figure  6.   To  first  order  in  this  system,  v.  =  x  cos  3.  - 
^  J  i  i 

y  sin  3-  where  x  and  y  are  the  unknown  target  coordinates  and 

9 .  -  d> .  =  (9.  -  3.)  -  (x  cos  3-  -  y  sin  3-)/r..   The  use  of  this 
11      11  i    2  11 

relation  implies  that  the  bearing  line  determined  by  3.,  the 

bearing  of  the  initial  estimate,  is  approximately  parallel  to  the 

true  bearing  line  determined  by  <$> .   .   The  use  of  this  relation  for 

all  stations  imposes  the  constraint  on  the  maximum  likelihood  bearing 

estimates  by  replacing  what  would  otherwise  have  been  estimates 

of  n  independent  bearings  4>  ,  through  <$>      by  estimates  of  two  independent 

quantities,  the  rectangular  coordinates  x  and  y.   It  also  implies 

that  the  initial  estimate's  range  from  a  station  is  approximately 

the  target's  range  from  a  station,  that  is,  that  the  initial 

estimate's  position  is  relatively  close  to  the  target's  position. 
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Since  u.  =  r.  (0.  -<$>.),    the  likelihood  of  observed  bearinqs 
1111  ^ 

8,,  0_,  ....  0   can  be  written  as: 

1  2  n 

n  1  1    n      2       2 

L(0         9    ,     .  ..,    6    )    =      n   = —  exp   -yl    uf/af 

1         l  n  1    /2tT   a.  ^    1      1      1 

i 

where    a.    =    r.e.    and   e.     is    the    standard   deviation   of    0..      The 
1111  i 

maximum    likelihood    estimates    for    <f>,    through    <b      are   determined 

1       3    n 

by  the  estimates  for  x  and  y  that  make  L(9,,  9_,  ....  9  )  a 

2  J  1    2         n 

maximum.   In  this  case,  making  L(9,,  9_,  .  ..,  9  )  a  maximum  is 

12         n 

n    2   2 
equivalent  to  making  1     (u./o.)    a  minimum.   So,  to  find  the  maximum 
^  ,11 

likelihood  estimates  x  and  y,  solve  the  following  two  equations 
for  x  and  y : 


3(^L)|     =  0    and    ii!!Lj!)_|     =  o 

X  =  X  X  =  X 

y  =  y  y  =  y 


With  the  constraint  given  by  u .  =  w.  -  x  cos  3.  +  y  sin  6.  where 

1       i  l                          l         J                i 

w.    =    r.     (9.    -    3.),    the    equations    can   be   written   as    follows: 
l  11  l 

n  2 

E     (w.    -    x   cos    3.    +   y    sin    3.)  (cos    3 . ) /o .    =    0 

,i                          l         2                i  li 


and 


11  2 

I     (w.    -    x    cos    3.    +    y    sin    3.)     (sin    3 . ) /a .    =    0. 

,i  l         J  l  li 
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In  terms  of  the  following  quantities: 


2       2  ? 

A  =  I(cos  $i)/oi    ,  b  =  Z(sin  6-  cos  6 • ) /o   , 


C  =  Z(sin2  3i)/o2  ,  D  =  I(w.  cos  6-)/a2  , 


E  =  Z (w.  sin  3- )/o2  , 
i      i    i  ' 


the  equations  become: 

Ax  -  By  =  D 

Bx  -  Cy  =  E 

The  solutions  are: 


2 


(1)  x  =  (BE  -  CD)/(B   -  AC) 
and 

(2)  y  =  (AE  -  BD)/(B2  -  AC) 

A  confidence  region  can  be  constructed  about  an  estimated 
position.   In  order  to  indicate  how  this  can  be  done,  a  proba- 
bility region  about  the  true  position  will  be  considered  first 

Both  x  and  y  are  values  of  random  variables.   If  a  new  set 
of  bearings  0, ,  0?,  . ..,  0   is  observed  (for  the  same  initial 
estimate  and  a  fixed  target),  in  general,  a  new  pair  of  values 
x  and  y  will  be  obtained. 
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If  X  and  Y  represent  these  random  variables,  from  (1)  and  (2)  , 

1       n       2 
(3)      x  =  — = I  (W./o.)  (B  sin  0.  -  C  cos  0.) 

(B  -AC)   1    X   1  1  1 


,       n      ~ 

(4)      Y  =  — i Z  (W./a7)(A  sin  0.  -  B  cos  0.) 

(B  -AC)   1    1   X  X  X 


where  W.  =  r.  (0.  -  0.  )  . 
l     11     l 

Since  X  and  Y  are  a  linear  combination  of  the  n  normal  random 
variables  W,,  W~ ,  .../  W  ,  or  equivalently  of  the  n  normal  random 
variables  Q,  ,  0-,  ...,  0    ,    they  have  a  joint  normal  distribution. 
Since  E(Wi)  =  ri(<j>i  -  8i)  ,   if  Bi  =  <j) .  for  i  =  1,  2,  ...,  n,  that 
is,  if  the  initial  estimate  of  the  target's  position  is  at  the 
target's  position,  E(W,)  =  0  for  i  =  1,  2,  ...,  n.   In  this  case 
E(X)  =  0  and  E (Y)  =  0  and  the  joint  normal  distribution  is  centered 
on  the  object's  position.   To  the  degree  of  the  approximations 
that  have  been  made,  this  is  also  true  if  the  initial  estimate  is 
not  at  the  target's  position. 

A  region  of  minimum  area  for  a  given  probability  of  contain- 
ment of  an  estimated  position  can  be  determined.   The  region  is 
bounded  by  an  ellipse  which  is  centered  on  the  object's  position 
and  whose  axes  lie  along  the  axes  of  an  x'y' -coordinate  system 
that  is  obtained  by  rotating  the  xy-coordinate  system  that  is 
centered  on  the  object's  position  through  an  angle  y.      In  this 
system,  a~,  ~,  is  0,  that  is  X'  and  Y1  are  independent  normal  ran- 
dom variables.   The  two  coordinate  systems  are  illustrated  in 
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Figure  7.   The  coordinates  of  a  point  in  the  two  systems  are 
related  by 


x'  =  x  cos  Y  -  y  sin  > 


y'  =  x  sin  y  +  y  cos  y 


These  relations  imply: 


2      2     2  2     2 

(5)    cu   =  o^    cos  y    -    2a        cos  y    sin  y    +   a      sin  y    , 

x1     x  xy  y 


/^\     2      2   .  2  2     2 

(6)    o^   =  a^  sin   y  +  2a^  cos  y  sin  y  +  o A  cos   y 

y'     x  Ay  y 


and 


2     2  2         2 

(7)    a^  ^   =  (a^  -  a  J  sin  y  cos  y  +  a_  (cos   y  -  sin   y 

x '  y '      x    y  xy 


where  y,  the  angle  of  rotation  of  the  coordinate  axes,  is  posi- 
tive in  the  clockwise  direction.   And  o~,~,  =  0  implies 

2a~~ 
tan  2y  =    -x^ 


2     2 

a      -    o 


With  the  initial  estimate  of  the  target's  position  at  the 

2 
target's  position  E(W.)  =  0  and  therefore   Var(W.)  =  a.   for 

i  =  1,  2,  ...,  n.   In  this  case  ,  from  (3)  and  (4) 

a^  =  — -^ *-   I  (1/a?)  (B  sin  B.  -  C  cos  6-)2, 

x    (B^-AC)Z   1 
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North 
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Figure  7.   An  elliptical  confidence  region  and  the  primed 
coordinate  system  in  which  the  covariance  a"i/vi 
is  zero.   The  center  of  the  ellipse  and  the 
origin  of  the  coordinate  systems  are  at  the 
target's  estimated  position.   The  estimated  bear- 
ing a  and  estimated  range  p  are  indicated  for  a 
reference  point.   The  major  axis  direction  is  6. 
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2        1     n      2  o 

o^    =   — j 7   Z    (l/a.)(A  sin  6.  -  B  cos  B-) 

y    (B  -AC)   1      1  x  x 


and 


a 


1 


xy    (B 


^ 9  I  d/a.)  (B  sin  6-  -  C  cos  6  • )  (A  sin  3-  -  B  cos  B-) 

-ACT  1  1  X  1  X 


Using  the  definition  for  A,  B  and  C,  the  above  become 


(8)  a2  =  — ^ 


x    (AC-B2) 


(9)  a2  =  — ^ 


(AC-B2) 


and 


B 


(10)  ~~         2 

V   ;  xy    (AC-B^) 


So  tan  2Y  =  2B/(A-C)  for  6-  =  <t>  •  ,  i  =  1,  2,  ...,    n. 

With  the  target's  position  known  and,  consequently,  4> .  known 

2    2 
for  1  =  1,  2,  ...,  n,  the  above  expressions  for  a-.  a~  ,o~-  and  y 

f  x  '   y  '  xy 

can  be  used,  since  the  initial  estimate  of  the  target's  position 

can  be  taken  as  the  target's  position. 

With  values  for  a~  ■  a-  ,  a--  and  y,    values  for  a-,  and  o~ 

x    y    xy  x        y 

can  be  found  by  using  equations  (5)  and  (6).   The  probability 

that  an  estimated  position  will  be  within  an  ellipse  of  semiaxes 

koA,  and  ko~ .  which  is  centered  on  the  target's  position  is 
x '        y ' 
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2 

1  -exp(-k  /2) .   This  result  can  be  found  by  integrating  the  bi- 

variate  normal  density  over  the  ellipse.   And  the  area  of  the 

ellipse  is   k  o~,o~,. 
c  x'  y ' 

Given  estimates  x  and  y  found  by  using  Equations  (1)  and  (2) , 

the  eLlipse  with  semi-axes  kaA,  and  ka~ ,  in  a  coordinate  system 

c-  x'        y  2 

that  is  centered  on  the  point  (x,y)  and  has  been  rotated  through 

2 

an  angle  y    is  a  1  -  exp(-k  /2)  confidence  region.   This  follows, 

since,  to  the  degree  of  the  approximations  involved,  the  bivar- 

iate  normal  distribution  of  X  and  Y  is  centered  on  the  target's 

2    2 

position.   The  confidence  ellipse  is  defined  if  a~  ,  a~  and  a~~ 

^  x  '  y      xy 

can  be  found,  that  is  if  the  elements  of  the  covariance  matrix 
can  be  found.   To  the  degree  of  the  approximations  involved, 
this  can  be  done  as  follows:   First,  assume  the  initial  estimate 

of  the  target's  position  is  at  the  target's  position.   Then, 

2    2  . 

values  for  a~   a~   a~^  and  y  can  be  determined  by  using  Equations 
x  i    y  i    xy  J  a   -a 

(8) ,  (9)  and  (10)   These  values  can  then  be  used  to  determine 

2     2 

a~,  ,  a-,  and  a~,~,  by  using  Equations  (5),  (6)  and  (7).   Now, 
x      y       x  y 

with  a  value  for  k,  a  confidence  region  can  be  constructed.   To 
the  degree  of  the  approximations  involved,  the  shape  of  the  con- 
fidence region  is  independent  of  both  the  target's  position  and 
of  the  initial  estimate  of  the  target's  position. 

For  the  case  where  bearings  are  taken  from  the  target  on  two 
or  more  stations,  0.  is  the  reciprocal  of  the  bearing  taken  from 
the  target. 

A  discussion  of  the  theory  of  bearings  only  position 
estimation  procedures  for  situations  similar  to  the  one  considered 
here  is  given  in  Reference  2. 
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The  following  equations  are  evaluated  in  the  program  to 
determine  the  coordinates  x*  and  y*  of  the  initial  estimate: 

x*  sin  (62-61)  =  [ p1  sin  (a, -6  ) ]  sin  6 


-  [  p2  sin  (ou-9 2)]  sin  6. 


and 


y*  sin  (62-ei)  =  [p   sin  (a,-9  )  cos  9„ 

-  [  P2  sin  ^a2~92^  COS  91  * 

Using  the  point  determined  by  two  lines  of  bearing  as 
the  initial  estimate  was  suggested  by  a  similar  procedure 
described  in  Reference  3. 
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APPENDIX  2:   Program  Descriptions 
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PEST  is  a  program  that  implements  a  position  estimation 
procedure  that  is  described  in  Appendix  1.   With  two  or  more 
observations  on  or  from  an  object  from  or  on  two  or  more 
stations,  the  program  will  generate  a  position  estimate  for  the 
object  with  an  associated  confidence  region.   After  program 
initiation,  the  program  user  is  presented  with  two  options:   (1) 
input  bearings  from  stations  on  an  object,  (2)  input  bearings  on 
stations  from  an  object.   To  choose  the  first  option,  press  T 
or  t.   To  choose  the  second  option,  press  S  or  s.   Bearing 
information  is  entered  in  the  following  order:   (1)  the  observed 
bearing  on  or  from  a  station,  (2)  the  bearing  of  the  station 
from  a  reference  location,  (3)  the  range  of  the  station  from  the 
reference  location  and  (4)  the  bearing  error.   The  reference 
location  can  be  a  station  location.   In  this  case,  the  station's 
bearing  and  range  from  the  reference  location  are  both  taken  to 
be  zero.   After  two  observed  bearing  observations  and  their 
associated  station  and  bearing  error  information  have  been  input, 
the  user  is  given  the  option  of  (1)  generating  an  estimate  or 
(2)  continuing  to  enter  bearing  observation  data.   A  position 
estimate  for  an  object  is  specified  in  terms  of  its  range  and 
bearing  from  the  reference  location.   An  associated  elliptical 
confidence  region  which  is  centered  on  the  estimate  is  specified 
in  terms  of  the  length  and  direction  of  its  major  axis,  the 
length  of  the  its  minor  axis  and  the  probability  (confidence) 
that  it  contains  the  object.   The  program  user  is  given  the 
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option  of  defining  the  confidence  region  in  terms  of  (1)  size 
or  (2)  containment  probability.   For  example,  the  size  of  an 
ellipse  with  a  containment  probability  of  .8647  is  2.   In  this 
case,  the  ellipse  is  referred  to  as  a  two  sigma  ellipse,  since 
the  semi-major  and  semi-minor  axes  of  the  ellipse  are  two  standard 
deviation  units  in  length.   To  choose  the  first  option,  press 
S  or  s.   To  choose  the  second  option,  press  P  or  p.   The  estimate 
can  be  recalled  by  pressing  E  or  e.   Additional  observations  can 
be  entered  by  pressing  C  or  c.   To  quit,  press  Q  or  q. 

The  position  estimation  procedure  requires  an  initial 
estimate  of  the  object's  position.   In  the  program,  the  initial 
estimate  is  at  the  intersection  of  the  two  bearing  lines  that 
correspond  to  the  first  two  bearings  that  are  input  to  the 
program.   Because  of  this,  the  first  two  bearing  inputs  should 
be  from  or  on  the  two  stations  that  are  estimated  to  provide  the 
best  initial  estimate.   Note,  if  bearing  errors  are  large 
relative  to  the  angular  separation  of  the  two  stations  as 
seen  from  the  object,  the  bearing  lines  from  the  two  stations 
may  not  intersect.   If  this  is  the  case,  the  reciprocal  bearing 
lines  will  intersect  and  a  gross  error  in  the  final  position 
estimate  could  occur. 

The  procedure  which  is  essentially  equivalent  to  one 
described  in  Reference  2  is  based  on  the  following  assumptions: 
(1)   Bearings  are  taken  on  or  from  an  object  of  unknown  position 
from  or  on  two  or  more  stations  of  known  position.   (2)  The 
distances  involved  are  such  that  the  object  and  the  stations 
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can  be  considered  to  be  located  on  a  plane  surface  (a  flat  earth) 
(3)   The  values  of  observed  bearings  on  or  from  stations  are 
determined  by  independent  normal  random  variables  each  with  a 
known  standard  deviation  (the  bearing  error)  and  a  mean  equal 
to  the  true  bearing  (zero  bias) .   (If  there  is  bias,  it  is  known 
and  removed.) 
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COMP  is  a  program  that  implements  a  procedure  for  combining 
a  group  of  position  estimates  for  an  object.   In  order  to  use  the 
procedure,  the  estimates  must  have  associated  elliptical 
confidence  regions.   When  used,  the  procedure  combines  the 
estimates  and  their  associated  confidence  regions  into  a  single 
composite  position  estimate  and  associated  elliptical  confidence 
region.   After  program  initiation,  the  program  user  is  prompted 
to  input  the  number  of  elliptical  areas  (confidence  regions)  to 
be  combined.   The  user  is  then  presented  with  two  options  for 
specifying  the  areas:   (1)  by  containment  probability  or  (2)  by 
size.   For  the  first  option,  the  areas  are  specified  by  the 
probability  (confidence)  that  they  contain  the  object.   This 
option  is  chosen  by  pressing  P  or  p.   For  the  second  option,  the 
areas  are  specified  in  terms  of  ellipse  size.   This  option  is 
chosen  by  pressing  K  or  k.   The  size  of  an  ellipse  with  a 
containment  probability  of  .8647  is  2 .   In  this  case,  the  ellipse 
is  referred  to  as  a  two  sigma  ellipse,  since  the  semi-major  and 
semi-minor  axes  of  the  ellipse  are  two  standard  deviation  units 
in  length.   The  orientation  of  an  ellipse  is  specified  in  terms 
of  the  direction  (the  angle  delta  in  the  program)  of  its  major 
axis.   The  position  estimates,  the  centers  of  the  elliptical 
areas,  are  specified  in  terms  of  latitude  and  longitude.   The 
user  is  given  the  option  of  determining  a  confidence  region  in 
terms  of  (1)  its  containment  probability  (confidence)  or  (2)  its 
size.   To  choose  the  first  option,  press  P  or  p.   To  choose  the 
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second  option,  press  K  or  k.   To  recall  the  position  estimate, 
press  E  or  e.   To  quit,  press  Q  or  q. 

The  procedure  which  is  described  in  detail  in  Section  III  is 
based  on  the  following  assumptions:   (1)   The  individual  estimates 
are  determined  by  independent  normally  distributed  random 
vectors  with  known  covariance  matrices  and  a  common  but  unknown 
mean  vector  whose  components  are  the  object's  rectangular 
coordinates.   (2)   The  object  is  located  in  the  plane  of  the 
coordinate  axes,  a  plane  tangent  to  a  spherical  earth.   In 
the  program,  the  coordinates  of  the  first  entry  are  the  coordi- 
nates of  the  point  of  tangency .   Because  of  other  uncertainties, 
this  flat  search  assumption  should  not  introduce  significant 
estimation  errors  for  the  distance  scale  for  which  the  program 
is  intended.   In  a  sense,  the  procedure  is  a  generalization  of 
the  procedure  described  in  Reference  4  which  is  limited  to 
circular  confidence  regions. 
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APPENDIX  3:   Program  Listings 
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PEST 

10  CLS:  SCREEN  0,0:  DIM  A<7) 

20  PRINT  "T=bearings  on  the  target.":  PRINT  "S=bearings  on  the 

stations. " 

25  A*=INKEY* 

30  IF  A*="T"  OR  A*="t"  THEN  JJ=0:  GOTO  60 

40  IF  A*="S"  OR  A$="s"  THEN  JJ=1:  GOTO  60  ELSE  GOTO  25 

60  CLS: INPUT  "observed  bearing" ;P:  P=P*ATN ( 1 > /45:  IF  JJ=1  THEN 

P=P+ 1 80 

70  INPUT  "station  bearing";Q:  Q=Q*ATN ( 1 ) /45: INPUT  "station 

range";R:  INPUT  "bearing  error";0:  0=0*ATN < 1 ) /45 

80  IF  1=2  THEN  GOTO  140 

90  1=1+1:  A(I-1)=P:  A(I+1)=Q:  A(I+3)=R:  A(I+5)=0:  IF  1=1  THEN 

GOTO  60 

100  X=A(4)*SIN<A(2)-A<0> ) :  Y=A (5) *SIN ( A (3) -A < 1 > > :  Z=SIN ( A ( 1 ) -A (O) 

) :  IF  Z=0  THEN  GOTO  330 

110  U=<X*SIN<  A(1)>-Y*SIN  (A(0)))/Z:  V=(X*COS  ( A ( 1 > ) -Y*COS (  A<0>> 

)  /Z 

120  FOR  M=0  TO  1 

130  P=A(M):  Q=A(M+2):  R=A(M+4):  0=A(M+6):  GOSUB  390:  NEXT  M:  GOTO 

150 

140  GOSUB  390 

150  CLS:  PRINT:  PRINT  "E=Est      C=Cont" 

160  A*=INKEY*:  IF  A$="E"  OR  A*="e"  THEN  GOTO  180 

170  IF  A*="C"  OR  A*="c"  THEN  GOTO  60  ELSE  GOTO  160 

180  CLS:  F=<B*B-A*C>:  IF  F=0  THEN  GOTO  380 

190  X1=U+(B*E-C*D)/F:  Y1=V+ ( A*E-B*D> /F:  GOSUB  440:  K=R1:  J=B1 

200  T=SGN(B)*ATN(1) :  IF  A=C  GOTO  220 

210  T=.5*ATN(2*B/ <A-C) > 

220  G=(C»C0S<T>*C0S(T>-2*B*C0S<T)*SIN<T>+A*SIN<T>*SIN<T>  >/-F: 

G=SOR(G) 

230  H=<C*SIN<T>*SIN<T>+2*B*C0S(T)*SIN(T>+A*C0S(T)*C0S(T> )/-F: 

H=SQR(H):  IF  H>=G  GOTO  250 

240  Z=H:  H=G:  G=Z:  T=T+2*ATN < 1 ) 

250  CLS:  PRINT:  PRINT  USING  "\         \####. ##" ; "bear i ng=" , J*45/ 

ATN  < 1 ) 

255  I1=K:  GOSUB  500:   PRINT  "range*" ; 1 1 s  PRINT 

260  PRINT  "S=Size  P=Prob  E=Est  C=Cont  Q=quit" 

265  A*=INKEY* 

270  IF  A*="S"  OR  A$="s"  THEN  CLS:  GOTO  320 

280  IF  A*="P"  OR  A*="p"  THEN  CLS:  GOTO  350 

290  IF  A$="E"  OR  A*="e"  THEN  CLS:  GOTO  250 

300  IF  A*="C"  OR  A$="c"  THEN  CLS:  GOTO  60 

310  IF  A$="Q"  OR  A*="q"  THEN  END  ELSE  GOTO  265 

320  INPUT  "size";S:  IF  S<=0  THEN  GOTO  320 

330  0=l-EXP(-S*S/2) 

340  PRINT  USING  "\  \##. ####"; "probabi 1 i ty=" ,0:  GOTO 

370 

350  INPUT  "probabi lity";0:  IF  0>=1  OR  0<=0  THEN  GOTO  350 

360  S=SQR(-2*L0G<1-Q) ) :  PRINT  USING  "\      \####. ##" ; "si ze=" ,S 

370  X=S*G:  I1=2*X:  GOSUB  500:PRINT  "major  axis=";Il 

371  N=T*45/ATN(1) :  IF  N<0  THEN  N=N+180 

372  PRINT  USING  "\  \####. ##"; "di recti on=" ,N 

375  Y=S*H: I1=2*Y:  GOSUB  500:  PRINT  "minor  axis=";Il:  I1=4*ATN(1)» 
X*Y:  GOSUB  500:  PRINT  "area=";Il:  GOTO  260 
380  PRINT  "no  solution":  END 
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90  X1=U-R*SIN<Q> :  Y 1=V-R*C0S (0) s  GO  3UB  440 

00  W=P -Bis  L=R1*0:  IF  L=0  THEN  GOTO  380 

10  G=C0S(B1) /L:  H=SIN(B1)/L:  IF  W>=4*ATN<1>  THEN  W=W-8*ATN ( 1 )  : 

OTO  43') 

20  IF  IaK=-4*ATN<1)  THEN  W=W  <-8*ATN  ( 1  ) 

30  W=W/0:  A=G*G+A:  B=G*H+B:  C=H*H+C:  D=W*G+D:  E=W*H+E:  RETURN 

40  R1=3QR<X1*X1+Y1*Y1) :  IF  R1=0  THEN  Bl=Os  RETURN 

50  IF  ^BS<X1/R1>=1  THEN  M1=SGN ( X 1 ) *ATN < 1 > *2  ELSE  M1=ATN ( X 1 /Rl / 

QR<1-X  L*X1/R1/R1> ) 

60  IF  .}BS(Y1/R1>=1  THEN  BI=2*ATN < 1 ) * < 1-SGN < Yl > >  ELSE  B1=2*ATN(1> 

ATN<Y1/R1/SQR(1-Y1*Y1/R1/R1 ) > 

70  IF  MKO  THEN  B1=8*ATN  ( 1 )  -Bl 

BO  RETURN 

00  11=100*11:  I1=INT(I1):  11=11/100 

10  RETURN 
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COMP 

5  CLS:  PI=4*ATN<1>:  INPUT  "Number  o-f  Elliptical  Areas" ;A0:  PRINT 

10  PRINT  "Area  De-fini  tions:  ":  PRINT 

11  PRINT  "By  Containment  Probability,  Press  P.":  PRINT 

12  PRINT  "By  Sigma  Size,  Press  K. " 

15  E*=INKEY* 

16  IF  E*="P"  OR  E*="p"  OR  E*="K"  OR  E*="k"  THEN  GOTO  17  ELSE  GOTO 
15 

17  CLS:  PRINT  "For  the  Latitude  and  Longitude  entry  -format,  press 
F.  Otherwise,  press  C. " 

18  K*«INKEY* 

19  IF  K*="C"  OR  K*="c"  THEN  GOTO  26 

20  IF  K*="F"  OR  K*="-f"  THEN  GOTO  23  ELSE  GOTO  IS 

23  CLS:  PRIN'I  "Latitude  and  Longitude  are  entered  in  degrees  and 
minutes  and  tenths  o-f  minutes  in  the  -form:   DDD-MM.MX   where  X  is 
N,S,W  or  E.   Leading  zeros  are  optional,  but  the   -   must  be 
included.   Press  C  to  continue. " 

24  C*=INKEY* 

25  IF  C*="C"  OR  C*="c"  THEN  GOTO  26  ELSE  GOTO  24 

26  CLS: FOR  0=1  TO  AO 

27  EE=0:  PRINT 

28  INPUT  "LAT";A$:  GOSUB  410:  IF  EE=1  THEN  GOTO  27 

29  Y=VA 

30  INPUT  "LONG"; A*:  GOSUB  410:  IF  EE=1  THEN  GOTO  27 

31  X=VA 

32  IF  0>1  THEN  GOTO  34 

33  Y0=Y:  X0=X:  C0=C0S ( Y0*PI / 180) 

34  Y=(Y-Y0)*60:  X=(X0-X) 

35  IF  X>270  THEN  X=X-360 

36  IF  X<-270  THEN  X=X+360 

37  X=X*60*C0 

38  INPUT  "Delta" ;S 

40  IF  S>=180  OR  S<0  THEN  GOTO  35 

45  INPUT  "Major  Axis"  ;P:  INPUT  "Minor  Axis";Q:  P=P/2:  Q=Q/2 

50  IF  E*="K"  OR  E$="k"  THEN  GOTO  60 

55  INPUT  "P";P3:  R=SQR (-2*L0G ( 1-P3) ) :  GOTO  70 

60  INPUT  "K";R 

70  IF  S>=135  THEN   S=S-180:T=P  :P=Q:  Q=T:  GOTO  110 

80  IF  S>=45  THEN  S=S-90:  GOTO  110 

90  IF  S>=0  THEN  T=P:  P=Q:  Q=T:  GOTO  110 

110  P=P/R:  0=0/ R 

120  S=S*PI/180 

1 30  U=P*P*C0S ( S ) *C0S ( S ) +Q*Q*S I N ( S ) *S I N ( S ) 

1 40  v=P*P*S I N ( S ) *S I N ( S ) +Q*Q*C0S ( S ) *C0S ( S ) 

150  W=(Q*Q-P*P>*SIN(S>*COS<S> 

160  Z=U*V-W*W:  A=A+V/Z:  B=B-W/Z:  C=C+U/Z:  D=D+V/Z*X-W/Z*Y 

170  E=E-W/Z*X  t-U/Z*Y:  F=F+V/Z*V/Z*U:  G=G+W/Z*W/Z*U 

180  H=H-V/Z*W/Z*U:  I=I+W/Z*W/Z*V:  J=J+U/Z*U/Z*V 

190  K=K+W/Z*U/Z*V:  L=L-V/Z*W/Z*W 

200  M=M+<WZ*U/Z+W/Z*W/Z)*W:  N=N-W/Z*U/Z*W 

210  NEXT  0 

215  CLS 

220  0=A*C-B*B:  X= (C*D-B*E) /0:  Y= <A*E-B*D) /0 

230  R= (C*C* (F+I+2*L) -2*C*B* (H+K+M) +B*B* (G+J+2*N) ) / (0*0) 
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j>40    S= (A*A* (G+J+2*N) -2*A*B* (H+K+M) +B*B* (F+I+2*L) ) / (0*0) 

[>50    T=(  <A*C+B*B)* (H+K+M) -C*B*<F+I+2*L)-B*A*(G+J+2*N)  )  /  (0*0) 

>55    Y=Y0+Y/60:     VA=Y:     GOSUB    470:     GOSUB    500 

>57    PRINT    "LAT    =    "  ;  L$ 

?5B    X=XO-X/ (60*C0) :     VA=X:     GOSUB    470:     GOSUB    540 

>59    PRINT    "LONG    =    " ; G* 

>70    A=SGN(T)*ATN(1) :     IF    S=R    GOTO    280 

>75    A=.5*ATN(2*T/(S-R) ) 

i?80    B=R*C0S(A)*C0S(A>-2*T*C0S(A)*SIN(A>+S*SIN(A>*SIN(A) 

?90    C=R*SIN(A>*SIN(A)+2*T*SIN(A>*C0S(A)+S#C0S<A)*C0S(A) :     A=A*180/ 

•I 

500    IF    B>=C    THEN    D=A+90:     GOTO    350 

505    E=B:     B=C:     C=E:     IF    A<0    THEN    D=A+180:     GOTO    350 

510    D=A:     GOTO    350 

!>20  PRINT:  INPUT  "P";P3:  K.=SQR  (-2*L0G  ( 1-P3)  )  :  Fl-li  GOTO  330 

525  PRINT:  INPUT  "K";K:  P3=1-EXP <-K*K/2> :  F1=0 

530  CLS:  PRINT  "Major  Ax i s  =  " ; 2*K*SQR <R) 

1535  PRINT  "Minor  Axis  =  " ;  2*K*SQR  (S) 

540  PRINT  "Delta  =  ";D:  IF  Fl=l  THEN  GOTO  345 

541  PRINT  "P  =  ";P3:  GOTO  350 
545  PRINT  "K  =  ";K 

K50  PRINT:  PRINT  "E  =  Est   P  =  P   K  =  K   Q  =  Quit" 

555  A*=INKEY* 

560  IF  A$="E"  OR  A*="e"  THEN  GOTO  390 

$65  IF  A*="P"  OR  A*="p"  THEN  GOTO  320 

570  IF  A*="K"  OR  A*- " k "  THEN  GOTO  325 

580  IF  A*="Q"  OR  A*="q"  THEN  GOTO  400  ELSE  GOTO  355 

590  CLS:  PRINT  "LAT  =  ";L*:  PRINT  "LONG  =  ";G*:  GOTO  350 

100  CLS:  END 

HO  VA=VAL(A*) 

U5  D*=MID*(A*,4,1)  ::  IF  D*="-"  THEN  GOTO  435 

120  D*=MID*(A*,3,1) :  IF  D*="-"  THEN  GOTO  430 

123  D*=MID*(A*,2,1> :  IF  D*="-"  THEN  GOTO  425 

t24  CLS:  PRINT  "Data  entry  error,  restart  entry  -for  the  area.": 

EE-1:  RETURN 

t25  M*=MID*(A*,3) :  GOTO  440 

V30  M$=MID$(A$,4) :  GOTO  440 

IZ5    M*=MID*(A*,5) :  GOTO  440 

140  VM=VAL(M$>:  VA=VA+VM/60 

[50    R*=RIGHT*(M*, 1) 

160  IF  R$="S"  OR  R*="s"  OR  R*="E"  OR  R*="e"  THEN  VA=-VA 

165  RETURN 

170  QQ=SGN(VA) :  VA=ABS(VA) 

173  IF  VA>180  THEN  VA=(VA-360)  ELSE  GOTO  480 

175  QQ=-QQ:  VA=ABS(VA) 

tBO    DD=F I  X ( VA ) :     FF=VA-DD :     MM=-F I  X ( 6000*FF ) / 1 00 

WO    RETURN 

500  IF  QQ=1  THEN  Q$="N":  GOTO  520 

10  IF  QQ=-1  THEN  Q*="S"  ELSE  Q*="  " 

20  L*=STR*(DD)+STR*(MM)+Q* 

30  RETURN 

40  IF  QQ=1  THEN  Q*="W":  GOTO  560 

50  IF  QQ=-1  THEN  Q*="E"  ELSE  Q*="  " 

60  G*=STR$(DD)+STR$(MM>+Q$ 

70  RETURN 
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